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A. J. Ahumada, Jr. and A. B. Watson 
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Introduction 


We have generated pictures which are samples of random 
contrast noise. The noise spectrum in a region of the picture 
varies with the distance from the region to the center. At 
each distance the noise spectrum is inversely proportional to 
the human spatial frequency contrast sensitivity function at 
the corresponding eccentricity. When an observer fixates the 
center of the picture from the appropriate distance the noise 
appears to have equal contrast at each spatial frequency and 
location. If contrast thresholds are set by internal noise, 
then this picture represents that noise. 


There is evidence that contrast sensitivity functions at 
different eccentricities differ only by a scale factor along the 
spatial frequency dimension. This means that such pictures can 
be generated by passing "white" noise through a space-invariant 
filter and then stretching the output by the appropriate 
space-varying magnification factor. The filter shape will be 
the inverse of the contrast sensitivity function at the fovea. 
The parameters for the filter and the magnification factor were 
taken from a model of spatial contrast vision devised by Watson 
(ref. 1) . 


The pictures summarize a noisy linear model of detection 
and discrimination of contrast signals by referring the internal 
noise of the model to the domain of the input picture. The 
detectability of an arbitrary pattern can be estimated by 
computing its ideal detectability in samples of the noise. 

Since unlimited inspection of contrast targets in noise leads to 
near ideal detection performance (Burgess et al., (ref. 2)), the 
absolute detectability of a low-contrast target in a particular 
retinal position for a brief duration can be estimated by adding 
a sample of noise to the target and trying to detect the target 
visually with no constraints on the viewing conditions. 



The Picture 


If one views figure 1 from a distance where it subtends four 
degrees of visual angle and fixates the center, the contrast 
variations in different regions of the picture should look roughly 
equal because they are all about the same level above threshold in 
log units. Both the contrast in spatial regions varying in 
eccentricity from the center and the contrast in different spatial 
frequency bands in each region have been made proportional to the 
threshold for detecting small contrast variations at that 
eccentricity and in that spatial frequency band. There is 
evidence that as actual contrast increases, there is a greater 
increase in apparent contrast for spatial frequencies with 
higher contrast thresholds (Georgeson and Sullivan, (ref. 3)), so 
the perceptual uniformity should be best close to threshold. 



Figure 1. A picture of the noise of the contrast detection 
system referred to its input. 


The perceptual uniformity was not the original goal in 
generating the picture. The original goal was to make a picture 
of noise equivalent to that visual system noise which limits our 
ability to detect spatial contrast signals. In engineering 
language, we wanted to make a picture representing the noise of 
the spatial contrast detection system referred to its input. 
Pelli (ref. 4) has discussed temporal and spatial charcter istics 
of this noise. Here we are concerned only with the spatial 
aspects of the noise which limits the detection of a brief 
(approximately one sec, with smooth onset and offset) 
presentation of a static, low-contrast signal. 
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Although the concept of referring the performance-limiting 
noise of a system to its input is common in engineering, it is 
not often encountered in perception. The advantage of reporting 
unidimensional absolute (or difference) thresholds in stimulus 
units was seen as soon as quantitative measurements of 
thresholds were attempted. This is equivalent to referring the 
noise of the observer back to the input. Consider these three 
statements about an observer's ability to detect brightness 
differences : 


(1) The observer has a brightness difference threshold of 
0.1 log units. 


(2) The observer has d' = l when the brightness difference 
is 0.1 log units. 


(3) The observer's noise referred to the input has a standard 
deviation of 0.1 log units. 


These three statements ascribe equivalent ability to the observer. 
Statements like (2) are currently the most popular formal 
expressions of thresholds. They result from models of the 
detection process in which the presentation of a stimulus value 
results in an internal random variable whose variability limits 
detection performance. Referring the noise to the input is just 
transforming the internal random variable back into the stimulus 
dimension. The appropriate transformation is the inverse of the 
transformation which computes the mean of the internal 
distribution from the value of the stimulus. 


In the case of multidimensional stimuli, the equivalent 
noise at the input will have a multivariate distribution in the 
dimensions of the stimulus. For the case of visual contrast 
detection for arbitrary targets, the stimulus can be regarded 
as a multidimensional stimulus with the contrast of each pixel 
being a dimension. The input equivalent noise will then be a 
multivariate distribution giving the probability that the pixels 
will jointly take on particular contrast values. The 
construction of the input equivalent noise for contrast 
detection is just the construction of a sample from that 
distribution. The theory behind the construction of the picture 
of the equivalent noise is outlined and the actual construction 
method is described. 
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The Theory 


The basis for the construction of the picture of equivalent 
noise is the signal-specif ied-exactly version of a model of visual 
contrast detection which has been shown to accurately predict a 
range of human contrast detection data by Watson (ref. 1) and 
Watson and Ahumada (ref. 7). In this model, contrast variations 
activate linear feature sensors whose receptive fields are similar 
to those of simple cortical cells. Noise is present at the output 
of these feature sensors. Both the filtering action of the 
feature sensors and the noise at their outputs restrict the 
ability of the model to detect low-contrast signals. The noisy 
output of the sensors is then combined linearly into a single 
number by the ideal decision mechanism, which then compares the 
number with a criterion do determine the response. Here we 
derive the equivalent input noise for this model. If this noise 
is added to any contrast signal and the sum is presented to an 
ideal observer, the performance level of the observer will be 
equivalent to that of the ideal observer looking at the noisy 
feature-sensor outputs. The equivalent noise thus reflects both 
the filtering and the noisy aspects of the model. 


Before showing the results for this model, we describe an 
apparently more general class of models for which the method is 
appropriate. The early stages of the visual system are organized 
in layers, so it would be convenient to model the visual system 
as a series of layers of processing. Figure 2 illustrates a 
general multilayer model for spatial contrast detection. A 
signal enters and is processed by a series of K layers. The 
spatial contrast signal is represented by a column vector S 
consisting of N contrast values, one for each pixel in the signal 
picture . 


S = ( s t ) , i = 1, N. 


( 1 ) 
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Figure 2. A multilayer model for contrast detection. 
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Pictures (spatial contrast signals) are usually represented as 
rectangular arrays, but for our purposes it is more convenient to 
place the rows end to end and regard a picture as a long row 
vector. In our actual computations, the picture resolution is 
256 by 256 pixels, so that N is 65,536. Each layer is represented 
as a noisy linear operator. For example, the first layer might 
represent the receptors, with the linear operator representing the 
optical point spread function, and the noise a combination of 
photon and transducer noise. In figure 2 the top boxes represent 
the linear operators, which are just matrices of coefficients. 


A k = ( a kij ) ' k = 1 , K; i = 1, N k _ 1? j = 1, N k (2) 

is the kth operator, where a ki j is the weighting that the jth 
output of the kth operator gives to its ith input and 


N 0 = N. (3) 


If the first layer represented the receptors, a^ij would reflect 

assumed to be Gaussian with a mean vector of zero and an 
arbitrary covariance matrix. 

^ k = ( t k i j ) , k = 1, K; i, j = 1, N k (4) 


is the kth covariance matrix, giving the covariance between the 
noise added to the ith and jth output variables at the kth stage. 
Returning to the receptor layer example, the off-diagonal 
elements of the noise covariance matrix would be zero to 
represent the independence of photon and receptor noise between 
two different receptors and the diagonal elements would have the 
variances of the outputs of the receptors. The output of the last 
stage goes to an ideal classifier, which, for the case of 
detecting the presence or absence of a single pattern, is a simple 
linear classifier. 


How reasonable are the above assumptions when the visual 
system is known to have nonlinear response functions and to have 
noise sources which are usually characterized as approximately 
Poisson, where the variance is a function of the signal strength? 
This is only a model for the detection of small contrast signals. 
Only the small signal response of the system has to be linear and 
the noise is assumed to be dominated by the noise generated by the 
background level itself. The ideal observer assumption can be 
regarded as a convenient way of describing the consequences of 
computing the noise of the system referred to its input. 
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Mathematically, the multiple layers can be represented by a 
single layer whose linear operator A is the product of the 
previous operators, 


A “ A^2 


a k . 


( 5 ) 


and whose noise distribution is computable from the noises and 
operators of the mulitple layer model* If the noise of the 
different layers is independent, the covariance matrix R for the 
equivalent single layer is given by 


R - R r+ A k \uT R k _ iA + ... + H^uT^ + ... + 
where 


H k = A k+l A k+2 * * • a K* ( 

Figure 3 illustrates this simplification. The two models are 
the same from the point of view of the ideal observer. In 
either case its inputs are normally distributed with mean SA and 
covariance matrix R. 



Figure 3. A single layer model for contrast detection. 


The performance level of the ideal observer as a function of 
the contrast signal S can be specified in terms of the hit rate 
and the false alarm rate by the measure d' given by 
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d'(S) = z (Prob{D I Signal=S} ) -z (Prob{D I Signal=0 } ) , 


( 8 ) 


where D is the observer's response that a signal was present 
and z is the functional inverse of the standard normal 
cumulative distribution function. The equation for the 
performance level of the ideal observer for the above models is 
given by 

d' (S) *■ S A R"” 1 A T S^ 1 . (9) 


The derivation of this equation can be found in standard texts 
such as Anderson (ref. 6, chapt. 12). 



Figure 4. A single layer model for contrast detection whose 
noise is "white". 


This model can be further simplified to a model with the 
same performance function, but with a simplified noise structure. 
In this new case, the noise added to each linear filter unit's 
output is an independent standardized z score with zero mean and 
unit variance. The operator in this new model, then, accounts for 
the effects of both the operators and the noises of the previous 
models. This version of the model is illustrated in figure 4, 
where it is called the "white" noise model because independent, 
identically distributed normal variates arranged in a rectangular 
array generate a picture in which all spatial frequencies have 
equal expected contrast energy. The equation for the linear 
operator C in this case in terms of the operator and noise of 
the previous model is given by 
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( 10 ) 


C = A B _1 , 


where B is a matrix of factor loadings of R, that is 

R ■ B t B. (11) 


Since this model is a special case of the previous one, the 
identity of performance can be verified by substituting C for A 
and the identity matrix I for R in equation 9. 



Figure 5. The input equivalent noise model for contrast 
detection . 


Figure 5 illustrates another model with the same 
performance function, but with no linear operator. In this 
version, the noise is in the same domain as the signal. Samples 
of the noise are pictures and the variances and correlations of 
the pixel values have all the information in the original 
operators and noises relevant to detection. The equation of the 
covariance matrix for the equivalent input noise in terms of the 
linear operator of the "white" noise model is simply 

R E = ( C C T ) -1 . (12) 

The identity of the performance function can be verified by 
substituting the identity matrix I for A and R fi for R in 

equation 9. The noise in this model is the noise of the 
previous models referred to the input picture domain. 
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Figure 6. The input equivalent model with "back-projected" 
"white" noise. 


The model developed by Watson (ref. 1) has the form of the 
"white" noise model. The values for the matrix C are available 
in the paper by Watson and Ahumada (ref. 5) . Figure 6 
illustrates a way of computing the equivalent noise when C is 
known. In this figure, the noise starts as independent z scores 
and then is "back projected" through an operator D which is the 
reverse (transpose) of the operator for the standard noise model 
multiplied by an inverse gain factor correction; that is. 


D * C T ( C C T ) _1 . 


(13) 


The factor after the transpose is called an inverse gain factor 
correction because in the case that C is an orthogonal 
transformation, it is a diagonal matrix containing the inverses 
of the squared lengths of the rows of C. In general, it also 
corrects for the correlations among the rows of C. 


The Calculation Method 


The initial idea for calculating the noise was just to use 
the transpose of C on "white" noise to obtain the picture. We 
realized that we had to correct for the gain factors, but we had 
not yet derived Equation 12 so we were not aware of the 
appropriate way to correct for the correlations. We generated 
low resolution pictures (64 x 64 pixels) using this method and 
described it in the ARVO abstract (Ahumada and Watson, (ref. 7)). 
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Since the model is not spatially homogeneous, we could not see 
any simple way to derive the gain factor correction matrix. 
However, the model is approximately a spatially homogeneous model 
preceded by a nonhomogeneous magnification function. That is, the 
linear dimensions of the linear feature detectors of the model are 
scaled by a factor 


s = 1 + ce , (14) 

where e is the eccentricity in degrees of the center of the 
response area from the center of the fovea and the constant c is 

estimated to be about 0.023 (mrad) - ^ (0.4 (deg. )”•*•). The linear 
dimensions of the detectors are thus doubled at an eccentricity 
of 44 mrad (2.5 deg.). The response areas of the detectors at 
the fovea range from about 0.7 mrad (0.04 deg.) to 91 mrad (5.2 
deg.) in octave steps. For all but the largest layers, the 
outputs of the detectors can be accurately approximated by first 
stretching the picture with a magnification function which is 
inversely proportional to that of equation 13, and then filtering 
with detectors whose response areas are independent of retinal 
eccentricity. This spatially homogeneous model can be represented 
by its response in the spatial frequency domain. Neglecting the 
presumeably small effects of the actual spatial position and 
orientation of the output features, the spatially homogeneous 
model can be approximated by a model which has uniform 
orientation and phase response, but with a spatial frequency 
amplitude response matched to that of the central fovea. 



Figure 7. "White" noise. The contrast value of each pixel 
is approximately normally distributed and independent of all 
other values. 


10 



Figures 7, 8, and 9 illustrate the actual calculations 
performed by the computer programs, whose listings appear in the 
appendix. All these programs were written in Fortran and were run 
under the RT-11 operating system on an LSI-11/23 processor. In 
addition to the program listings, the appendix also contains the 
operating system commands and program parameters used to create 
the actual files of data displayed in the figures. Like figure 1, 
figures 7 and 9 were photographed from the monitor of our raster 
graphics display. 


Figure 7 shows the "white" noise with which we started. It 
was generated by the program FNOIS, which computes an 
approximately normally distributed value for each pixel by 
adding 12 uniformly distributed pseudorandom numbers. We then 
transformed this picture into the spatial frequency domain using 
a floating point two-dimensional FFT program, FFFT, Next, the 
noise in the transform domain was multiplied by the inverse of an 
approximation to the spatial spectral sensitivity of the fovea. 
The sensitivity function used is given by 


s ( f ) = exp(-(f/14.) 2 ) - 0.9 exp(-(f/1.0) 2 ) , (15) 


where f is the spatial frequency in cycles per degree. This 
function is plotted in figure 8A. The amplitude of the noise in 
the frequency domain was multiplied by the factor a(f) given by 


a (f ) = s(32)/s(f) , f < 32 

= 0, f > or = 32. (16) 


This function is plotted in figure 8B. Since the 256 pixels 
across a picture are assumed to span 4 degrees of visual angle, 
frequencies in either the horizontal or vertical directions are 
limited to 32 cycles per degree. Higher frequencies do occur in 
diagonal directions, since the spatial frequency is given by 


f " ( f x\ u2 + f y \u2 j 1/2 ' 


(17) 


in terms of the frequencies in the x and y directions alone. The 
amplitude of the noise at these frequencies was limited so that 
essentially invisible noise power at these frequencies would not 
use up the small dynamic range of the output, which was limited 
to 8 bits. The frequency domain filtering was done by the program 
FDOGM. 
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SPATIAL FREQUENCY, cycles/deg 

Figure 8. (A) The model's contrast sensitivity at zero 

eccentricity as a function of spatial frequency. (B) The 
amplitude factor of the filter used to make the homogeneous 
filtered noise from the white noise. 


Figure 9 shows the result of transforming the filtered noise 
back into the spatial domain. As figure 8B indicates, this is 
essentially high-pass noise. We expected it to be more like 
band-reject noise because the spatial frequency sensitivity 
function is usually plotted in logarithmic frequency co-ordinates 
and the sensitivity function looks reasonably symmetric in these 
co-ordinates. 


Finally, the noise of figure 1 was obtained from the noise 
shown in figure 9 by stretching the noise according to the 
inverse of the cortical magnification factor. Using the 
subscript 1 for the spatially homogeneous input picture and the 
subscript 2 for the stretched output picture, the intensity I for 
the output is just copied from a corresponding intensity in the 
input, that is 


I 2 < x 2 'Y 2 ) = I i( x i»yi) • (IB) 
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r, 


Figure 9. The spatially homogeneous filtered noise. 

The input co-ordinates in degrees of visual angle corresponding 
to an output were found by 

= m x 2 and y-^ = m y 2 , (19) 

where the inverse expansion factor m is given by 



m = In ( 1 + cr 2 )/ C r 2 , 


( 20 ) 


where r 2 is the length of (x 2 ,y 2 ), t ^ iat i s 


r 2 - ( x 2 \u 2 + y 2 \ u 2 ) i/ 2 . 

The constant c is the parameter in the inverse local cortical 
magnification function in Watson's model (ref. 1) given by 


( 21 ) 


M (e) =1/ (1+ce) , 


( 22 ) 


which is just the multiplicative inverse of equation 13. 
Equation 20 is just the integral of this local magnification 
function along the radius to the point. The inverse of the 
constant c is the eccentricity in degrees at which the inverse 
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local cortical magnification function has dropped to a value of 
one half. We set this value to 44 mrad (2.5 deg.), which 

— 1 —1 

corresponds to a value of 0.023 (mrad) x (0.4 (deg.) A ) for c. 

There is an additional complexity in the program BCRMG which 
stretched figure 9 to obtain figure 1. The co-ordinates for a 
point in the input corresponding to a point in the output do not 
usually fall exactly on a point in the input. The output in the 
general case has to be interpolated from the surrounding points. 
In this case we used a four by four interpolation function which 
is the product of the sine function in each direction, with the 
amplitude of the interpolation normalized to correct for the fact 
that the sine function was only computed at the four closest 
values in each direction. 


Summary 


To summarize, the noise of figure 1 is perceptually 
homogeneous in the spatial and spatial frequency domains. Except 
for a contrast gain factor, it represents in a single picture the 
spatial contrast detection ability of the visual system. It is a 
picture of the noise of the visual contrast detection system in a 
useful sense. Once any signal is added to the noise the result 
can be inspected at any magnification, gain, or for any length of 
time to see whether the signal is "visible". 
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APPENDIX 


This appendix contains the Fortran IV source code listings of 
the programs used to create the figures, preceded by a list of the 
required RT-11 system commands and program control parameters. 


1 RT-11 system commands and program parameters 
! Device m: should be fastest device available 
I (we used a ramdisk device driver), 
create m: image. flo [1024] 
r fnois 

create image . byt [ 128] 
r fb 

28.0,0.0 

copy image. byt noise. byt 
I figure 7 in file noise. byt 
r ffft 
0 

r fdogm 
r ffft 
-1 

r f b 

100 . 0 , 0.0 

copy image. byt noidog.byt 
1 figure 9 in file noidog.byt 
copy image. byt m: input. byt 
r bcrmg 

1 figure 1 in file image. byt 


program fnois 

c This program fills the file M: IMAGE. FLO with 256 times 256 
c floating point numbers which are approximately independent, 
c normally distributed values with a mean of 0.0 and a 
c standard deviation of 1.0. These numbers are followed by an 
c equal number of zeros in anticipation of a complex fft. 
c 

real z, row(256) 
c row, column dimensions 
n= 256 

open( unit=l, name= 'MrlMAGE.FLO' , type='OLD', 

& access= 'DIRECT', recordsize= n) 
c random seeds 
iran= 0 
jran= 0 
c do rows 

do 6 j=l, n 
c do columns 
do 7 i= 1, n 
z= 0 . 

do 10 k= 1, 12 
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10 z= z + r an ( iran , j r an) 
row ( i ) = z - 6.0 

7 continue 

jl= j 

write(l'jl) (row(i), i= 1, n) 
type *, j, row(l), row(n) 

6 continue 
c do zeros 

do 8 i= 1, n 

8 row ( i) = 0.0 
do 9 j= 1, n 
jl= j + n 

9 write(l'jl) (row(i), i= 1, n) 
stop 

end 


program fb 

c Converts floating input of n records of length n, 
c to byte output (-128 to 127) using scale factors input 
c from the terminal. Scale factors are obtained from 
c program frang. 
c 

real c(256) 
byte a(256) 
integer ia 
n= 256 

open ( unit= 1, name= ' M: IMAGE. FLO' , type= 'OLD', 

& access= 'DIRECT', recordsize= n) 
nd4- n/4 

open ( unit= 2, name= 'IMAGE. BYT' , type= 'OLD' , 

& access= 'DIRECT', recordsize= nd4) 
type*, ' scale factor, offset ?' 
accept*, scale, offset 
do 7 j= 1, n 

jl= j 

read(l'jl) (c(i), i=l, n) 
do 8 i= 1, n 
fa = scale*c ( i) +of f set 
if ( fa .It. 0. ) fa = fa - 1. 
a ( i) = f a 
8 continue 

write(2'jl) (a(i), i=l, n) 
type* , ' row ' , j 1 
7 continue 
stop 
end 
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program frang 

c Finds the maximum and minimum values in a complex image, 
c 

real c (256) , d (256) 
n= 256 

open ( unit=l, name= 'M: IMAGE. FLO ' , type= 'OLD', 

& access= 'DIRECT', recordsize= n) 

cmax=-10000000 . 

cmin= 10000000. 

dmax=cmax 

dmin=cmin 

do 7 j= 1, n 

jl= j 

jml= j 1+256 

read (1 ' jl) (c (i) , i= 1, n) 
read (1 ' jml) (d ( i) , i= 1, n) 
do 8 i= 1, n 
cmax=amaxl (cmax,c (i) ) 
dmax=amaxl (dmax ,d ( i) ) 
cmin=aminl (cmin,c (i) ) 
dmin=aminl (dmin,d (i) ) 

8 continue 

type* , ' row ' , jl 
7 continue 

type*, ' cmax, cmin, dmax, dmin=' ,cmax,cmin, dmax, dmin 

stop 

end 


program ffft 

c Two-dimensional complex floating point fft. 
c Inverse flag is 0 for direct, nonzero for inverse, 
c Calls one-dimensional subroutine: 
c fft( n, c, d, costab, invers) 

c Output is written over input file, M: IMAGE .FLO 
c 

real c(256), d(256), costab(65) 
real cbuf(8,256), dbuf (8,256) 
n= 256 
n2= n*2 
lbuf = 8 

open( unit=l, name= ' M: IMAGE . FLO ' , type='OLD', 

& access= 'DIRECT', recordsize= n) 
pi = 3.1415927 
twopi =2. * pi 

type*, ' non-zero entry for inverse ? ' 
accept 100, invers 
100 format(ilO) 

if( invers .ne. 0 ) inverse= -1 
c fill cosine table 
do 10 i= 1, n/4 +1 
10 costab(i) = cos ( twopi* ( i-1 ) /n ) 
c type*, ( costab(i), i= 1, n/4 +1) 
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c do row transforms first 
do 1 j= 1, n 

jl= j 

jml= jl+ n 

read(l'jl) (c(i), i= 1, n) 
read(l'jml) (d(i), i=l , n) 
call fft( n, c, d, costab, invers) 

j 2= j 

jm2= j2+ n 

write(l'j2) (c(i), i = 1, n) 
write(l'jm2) (d(i), i= 1, n) 
type *, ' row j 

I continue 
close (unit=l) 

c do columns 
nbuf= n/lbuf 

open ( unit=l, name= 'M: IMAGE. FLO' , type= 

& access= 'DIRECT', recordsize= lbuf) 

nrec= n* (n/lbuf) 

do 6 il = 1, nbuf 

ibuf = il 

do 9 j= 1, n 

jl = ibuf 

ibuf = ibuf + nbuf 
jml = jl + nrec 

read(l'jl) (cbuf(i,j), i= 1, lbuf) 
read(l'jml) (dbuf(i,j), i=l, lbuf) 

9 continue 

do 11 i2 = 1, lbuf 
do 7 j= 1, n 
c (j ) = cbuf ( i 2 , j ) 
d ( j) = dbuf ( i 2 , j ) 

7 continue 

call fft( n, c, d, costab, invers) 
do 8 j = 1, n 
cbuf (i2,j)= c ( j ) 
dbuf ( i 2 , j)= d ( j ) 

8 continue 

icol= i2+ lbuf* ( il-1) 
type *, ' col ', icol 

II continue 
ibuf= il 

do 12 j= 1, n 
jl= ibuf 

ibuf= ibuf + nbuf 
jml= jl + nrec 

write(l'jl) (cbuf(i,j), i= 1, lbuf) 
write(l'jml) (dbuf(i,j), i=l, lbuf) 

12 continue 
6 continue 
stop 
end 




'OLD' , 
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program fdogm 

c Multiplies the amplitudes of a complex frequency domain 
c image by an amplitude response having the form of the 
c inverse of the difference of two Gaussians. 

c The 256 pixels are assumed to subtend 4 degrees of visual angle, 
c Spatial frequencies above 32 cycles per degree are given 
c an amplitude of zero, 
c 

real c, d 

real crow(256), drow(256) 
real gl( 256) , g2( 256) 
n= 256 

open ( unit= 1, name= ' M : IMAGE . FLO 1 , type= 'OLD', 

& access= 'DIRECT' , recordsize= n) 
nd2= n/2 
nd2pl= nd2+l 

c Difference of Gaussian parameters: ampl , spreadl , amp2 , spread2 
c are respectively, al, bl, a2, b2. 
al=l. 
a2= . 9 

bl= 10. *4. 
b2= 0.7*4. 
nx= 3.5*bl 
nx=minO ( nd2pl , nx) 
c= -0 . 5/ (bl*bl) 
do 10 i= 1, nx 
x = i-l 

10 gl(i)= al*exp (c*x*x) 
if(nx .ge. nd2pl) goto 14 
do 11 i=nx+l,nd2pl 

11 gl ( i) = 0. 

14 continue 
nx= 3.5*b2 
nx=min0 ( nd2pl, nx) 
c= -0.5/(b2*b2) 

do 12 i= 1, nx 
x=i-l 

12 g2(i)= a2*exp (c*x*x) 
if(nx .ge. nd2pl) goto 15 
do 13 i= nx+l,nd2pl 

13 g2 ( i) =0 . 

15 continue 

nx= minO ( 32, nd2pl) 
dogmin= al*al*exp (- . 5* (32/10. ) **2) 
type* , dogmin 
c do rows 

do 6 j=l, n 

jl= j 

read(l'jl) (crow(i), i=l, n) 
jl= j+ n 

read(l'jl) (drow(i), i=l, n) 
c do columns 
do 7 i= 1, n 
ix=i 

if ( ix .gt. nd2pl) ix= n-ix+2 
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iy= j 

if ( iy .gt. nd2pl) iy= n-iy+2 

factor= dogmin/( gl ( ix) *gl ( iy) -g2 ( ix) *g2 ( iy) ) 
if (factor .gt.l . ) factor= 0. 

8 crow ( i) =f actor*crow ( i) 
drow(i) = factor *drow ( i) 

7 continue 

jl= j 

write(l'jl) (crow(i), i= 1, n) 
jl= j + n 

write(l'jl) (drow(i) , i= 1, n) 
type* , 1 row 1 , j 
6 continue 
stop 
end 


program bcrmg 

c This program generates a picture in IMAGE. BYT 
c by sampling and interpolation from M: INPUT. BYT 
c according to the linear approximation to the cortical 
c magnification factor. 

c A sine function is used for the interpolation, with a gain 
c correction factor to compensate for amplitude variations 
c caused by truncation of the sine function, 
c 

byte row2 (256 ) , bufl(256,4) 

real w(201) 

real cor (100 ) 

n= 256 

nd4= n/4 

open ( unit= 1, name= ' M : INPUT. BYT ' , recordsize= nd4, 

& access= 'DIRECT', type= 'OLD') 

open ( unit= 2, name= 'IMAGE. BYT', recordsize= nd4, 

& access= 'DIRECT', type= 'OLD') 
pi = 3.1415927 
twopi =2. * pi 
c Sine weighting function, 
do 200 i= 2, 201 
il= i-1 

200 w(i)=sin(pi*il/100.)/(pi*il/100.) 

w(l)=l. 

c Amplitude correction, 
do 300 i= 1, 100 

300 cor ( i) =1 . 0/ (w ( i) +w ( i+100 ) +w (102-i) +w ( 202-i) ) 

c viewing distance in pixels so 1 pixel is 1/64 degree 
vdist=57*64 

c coefficient of magnification factor in pixels 
c= 0.4 * 57. / vdist 
ci= 1. / c 

c lower and upper ranges of output pixel ranges 
nl2 = 1 
nu2 = n 

c input and output magnification center coordinates 
no2 = ( nu2 + nl2 ) / 2 
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origin = no2 
jold= -1 

c For each output row 
do 10 j2 = nl2 , nu2 

c For each output column 
do 1 i2 = nl2 , nu2 
di = i2 - no2 
d j = j 2 — no2 

r2 = sqrt ( di * di + dj * dj ) 
if ( r2 .eq. 0. ) goto 2 
rl = ci * alog ( 1. + c * r2 ) 
ratio = rl / r2 
goto 3 

2 ratio = 0. 

3 

yl = 
il = 

jl = 

idx= 
idy= 
wx2= 
wy2 = 
wxl = 
wyl = 
wx3 = 
wy3 = 
wx4 = 
wy4 = 
ilpl 
jlpl 
ilml 
jlml = jl - 1 
ilp2 = il + 2 
j lp2 = jl + 2 
jrecl= jlml 

if (jlml .eq. jold) goto 9 
jold= jlml 
do 8 j= 1, 4 

read (1 ' jrecl) ( bufl(i,j), i=l,n ) 

8 jrecl= jrecl+1 

9 continue 


all 

= 

buf 1 ( 

ilml , 

1) 

al2 

= 

buf 1 ( 

ilml , 

2) 

al3 

= 

buf 1 ( 

ilml , 

3) 

al4 

= 

buf 1 ( 

ilml , 

4) 

a21 

= 

buf 1 ( 

il. 

1) 

a22 

= 

buf 1 ( 

il, 

2) 

a23 


buf 1 ( 

il, 

3) 

a24 

= 

buf 1 ( 

il, 

4) 

a31 

= 

buf 1 ( 

ilpl. 

1) 

a32 


buf 1 ( 

ilpl. 

2) 

a33 

= 

buf 1 ( 

ilpl , 

3) 

a34 

= 

buf 1 ( 

ilpl, 

4) 

a41 

= 

buf 1 ( 

ilp2 , 

1) 


xl = ratio * di + origin 
ratio * dj + origin 
xl 

yi 

int ( ( xl-il) *100 . ) +1 
int ( (yl- j 1) *100. ) +1 
w ( idx) 
w(idy) 
w (100 + idx) 
w (100+idy ) 
w (102-idx) 
w (102-idy) 
w (202-idx) 
w (202-idy) 

= il + 1 
= jl + 1 
= il - 1 
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a42 = buf 1 ( ilp2 , 2) 
a43 = buf 1 ( ilp2 , 3) 
a44 = bufl( ilp2, 4) 

ave = wxl* (wyl*all+wy2*al2+wy3*al3+wy4*al4) 
& +wx2* (wyl*a21+wy2*a22+wy3*a23+wy4*a24 ) 

St +wx3* (wyl*a31+wy2*a32+wy3*a33+wy4*a34) 

St +wx4* (wyl*a41+wy2*a42+wy3*a43+wy4*a44 ) 

row2(i2) = ave*cor ( idx) *cor ( idy) 

1 continue 
jrec= j 2 

write (2 ' jrec) (row2(i), i= 1, n) 
type* , ' row 1 , jrec 
10 continue 
stop 
end 
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